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Abstract.  The  effective  speed  of  sound  is  calculated  for  a  medium 
containing  immovable  rigid  spheres  arranged  in  a  simple  cubic  lattice. 
Long  waves  propagating  along  a  lattice  axis  are  treated.  The  wave 
equation  for  the  pressure  is  reduced  to  an  ordinary  differential  equation 
to  which  Floquet  theory  is  applied.  Both  perturbation  and  numerical 
methods  are  use  to  find  the  effective  speed  as  a  function  of  frequency, 
and  to  locate  the  pass  and  stop  bands. 
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1.  Introduction 


The  effective  speed  of  sound  in  a  medium  containing  particles  is  different  from  that  in 
the  ambient  medium.  We  calculate  it  for  the  idealized  case  of  immovable  rigid  spheres 
with  centers  arranged  in  a  simple  cubic  lattice.  We  assume  that  the  wavelength  in  the 
ambient  medium  is  large  compared  to  the  distance  between  particles.  Then  the  wave 
equation  for  che  pressure  can  be  reduced  to  an  ordinary  differential  equation  with  periodic 
coefficients.  Floquet  theory  shows  that  there  are  pass  bands  and  stop  bands  along  the 
frequency  axis.  We  determine  them,  and  the  sound  speed  in  the  pass  bands,  by  both 
analytical  and  numerical  means. 

2.  Formulation 

Let  the  x-axis  be  a  lattice  axis  with  L  the  spacing  between  sphere  centers  along  it,  and 
let  R  be  the  radius  of  each  sphere.  A  wave  of  frequency  u*  propagates  along  this  axis. 
By  the  periodicity  of  the  lattice,  the  wave  is  symmetric  about  the  planes  y  =  ±1/2  and 
z  —  ±L/2.  Therefore  it  suffices  to  determine  the  wave  within  the  region  bounded  by  these 
four  planes,  which  may  be  thought  of  as  a  rigid-walled  waveguide  of  square  cross  section 
with  spheres  placed  periodically  along  its  axis.  When  ujL/c  is  small,  where  c  is  the  sound 
speed  in  the  ambient  medium,  the  wavelength  in  this  medium  is  large  compared  to  I. 
Then  the  pressure  p  is  practically  constant  over  the  cross-section  of  the  waveguide,  so  we 
write  it  as  p(x). 

Under  these  conditions,  p(x)  satisfies  the  long  wave  equation  [Ll] 
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Pzt  +  ~jPt  +  k2p  =  0,  k-ui/c.  (2.1) 

Here  S(x)  is  the  unobstructed  cross-sectional  area  of  the  guide,  given  by 


S(x  +  L)  =  S(x). 


In  view  of  (2.2),  the  ordinary  differential  equation  (2.1,  has  a  periodic  coefficient.  Ey 
Floquet's  theorem  it  has  a  solution  which  satisfies  the  condition 

p(x  +  L)  =  t'KLp{x)  (2.3'' 

for  some  constant  K(k)  which  depends  upon  k.  The  complex  conjugate  pm  of  p  is  another 
solution  which  satisfies  (2.1)  with  iK  replaced  by  —iK’.  When  K  is  real,  (  2.3)  represents 
a  wave  propagating  in  the  direction  of  increasing  x  with  wavenumber  K .  Therefore  its 
phase  velocity  is  C  =  u/K  and  its  group  velocity  is  Cg  =  (dK/du:)-1 .  In  terms  of  k  these 
relations  become 

c=k_  £1_^diL=(d_Kyl 

c  K'  c  dK  \dk  )  ' 

When  K  is  not  real  p  increases  or  decreases  exponentially  with  x  and  the  wave  is  called 
evanescent  or  nonpropagating.  In  this  case  we  say  that  u  and  k  lie  in  a  stop  band,  while 
when  K  is  real  they  lie  in  a  pass  band. 

Thus  our  task  is  to  find  a  solution  p  of  (2.1)  and  a  constant  K(k)  such  that  (2.3) 
holds.  Instead  we  shall  specify  K  and  seek  k(K)  and  p  satisfying  (2.1)  and  (2.3)  because 
this  is  a  standard  eigenvalue  problem  with  eigenvalue  k2. 

It  will  be  convenient  to  introduce  dimensionless  variables  with  L  as  the  unit  of  length. 
Thus  we  set  x'  =  x/Z,,  R‘  —  R/L,  k1  =  kL,  K'  —  KL  and  S'  =  S/L2.  From  now  on  we 
will  use  these  variables,  omitting  the  prime. 

3.  Perturbation  expansion 

The  coefficient  of  pz  in  (2.1)  is  Sz/S  =  h rx/S  for  |x|  <  R  and  zero  for  R  <  lx|  <  1/2. 
Its  maximum  absolute  value  is  tR  w'hen  R  <  (27r)-$,  so  it  is  small  when  R  is  small 
Therefore  we  shall  solve  the  problem  by  a  perturbation  expansion,  treating  Sz/S  as  small 
The  results  will  be  valid  for  R  small,  although  we  shall  see  that  they  are  useful  even  for 
R  =  1/2,  when  adjacent  spheres  touch  one  another. 
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A  O  -O  HA 


To  proceed  we  write 


p(x)  =  Pol1)  +  Pi(x)  +  pi(x)  +  •  •  ■ 
k7  =  +  (^2)l  +  (^2)2  +  ■  •  •  ■ 


(3.1) 


A  term  with  subscript  j  is  of  jth  order  in  5r/5.  Now  we  substitute  (3.1)  into  (2.1)  and 
equate  terms  of  each  order  to  obtain 


Po  +  koPo  =  0* 

(3.2) 

p"  +  ^'oPi  =  ~yPo  -  (fc2)iPo< 

(3.3) 

p'i  +  kin  =  -jp'i  ~  (k2)m  -  2(iJ),p,. 

(3.4) 

We  also  find  that  each  p}  satisfies  (2.3). 

By  writing  po  =  eikaX  and  then  using  (2.3)  to  determine  k0  we  find  an  infinite  number 
of  solutions: 

po  =  e'koX,  kQ  —  K  +  2irm,  m  =  0,  ±1, -  (3.5) 

When  K  is  real  then  Jt0  is  real.  Before  solving  (3.3)  we  multiply  it  by  pi  -  e~'k°x  and 
integrate  the  resulting  equation  from  -1/2  to  1/2.  After  integrating  by  parts  and  using 
(2.3)  we  find  that  the  integral  of  the  left  side  vanishes.  Thus  the  integral  of  pi  times  the 
right  side  vanishes,  so  we  can  solve  it  for  (fc2)i  and  obtain 

(*’),  = -i/ho  /‘/2  =  0.  (3.6) 

J  —  1/2  --H1) 

Next  we  proceed  in  the  same  way  with  (3.4)  and  again  the  integral  of  p0  times  the 
left  side  vanishes.  Then  the  integral  of  pi  times  the  right  side  yields 

(*2)2  =  -  f  e~ik°ZJ(Z)pi(x)dx‘ 

J- 1/2  *-Hz) 

Now  we  solve  (3.3)  by  means  of  a  modified  Green’s  function  G(x,y),  to  be  determined 
below,  and  we  get 

p,(x)  =  -ih  fU1  (3.8) 

J- 1/2 
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Then  we  use  (3.8)  in  (3.7)  and  find 


(k\  =  ,k t.  /  f  3.9 1 

J-MiJ-\n  S(x)  5(y) 

The  results  (3.5),  (3.6)  and  (3.9)  yield  the  first  three  terms  in  the  expansion  of  k2 ,  while 
(3.5)  and  (3.S)  yield  two  terms  in  the  expansion  of  p(x). 

To  complete  the  calculation  we  must  determine  G  which  is  defined  by 


GtI  +  k\G  =  S(x  -  y)  -  eiko(t-v),  -1/2  <  x  <  1/2 


G(l/2,y)  =  e,KG(-l/2,y), 
Gx(l/2,y)  =  «■*£,(  — 1/2, y). 


The  solution  of  (3.10)  is  not  unique  when  ko  satisfies  (3.5).  One  solution  is 


/3. 10) 


eiio(i-y) 


x£ik0(z-y) 


2iko  2ikQ(l^e2'K)  2ik0  ’ 

£ik0(y-z)  xeik0(z~y) 


y<x<  1/2, 


(3.11) 


2ifc0(l  -  e2'K)  2ik0 


—  1/2  <  x  <  y. 


By  using  (3.11)  for  G  in  (3.8)  and  (3.9)  we  can  calculate  pi(x)  and  (fc2)2- 

Instead  of  (3.11)  we  shall  represent  G,  in  terms  of  the  eigenfunctions  (3.5): 

Gn(x,y)  =  ]T  ei(K+2*rn)(z-v)  ^  +  2^  -  +  2^  (312) 

m^n 

Here  we  have  set  k0  =  K  4-  27rn  and  denoted  G  by  Gn  We  shall  also  denote  p:  by  pni  and 
(k2h  by  (k2)n 2.  Then  (3.8)  yields 

pnl(x)  =  -l^e*^+2,rm)lVnm  (k  +  2™)  -  (k +2*m^j  (3.13) 


The  matrix  elements  Vnm  are  defined  by 


■  i: 


1/2  S'(y)  _,2ir(  n  — m)y  . 


(3.14) 


Similarly  (3.9)  yields 


(k2)n2  =  -(A  +  2?rn)  £  (A  Hr  27rm)|Vnm|2[(A  +  2?r n)2  -  (A  +  2rrm)2]-1.  (3.15) 

m#n 

Let  us  now  set  n  =  0  and  simplify  (3.15)  in  the  long  wavelength  case  when  A  is  small. 


We  obtain 


€  £iv“' 


m^O 


^-1 
m  — 


A 


2? rm2 


+  0(A2) 


The  first  sum  vanishes,  so  we  get 


K2 


(*2)o2  =  ~2^  E  +  °(K2y 


(3.16) 


(3.17) 


m>  1 


From  (3.5)  see  that  fcoo  =  A,  while  (3.6)  yields  (fc2)oi  =  0.  We  now  use  these  values  and 
(3. 1*7)  in  (3.1)  to  get  k2 .  Then  we  take  the  square  root,  divide  by  A  and  use  (2.4)  to  write 


C(R ) 


k_ 

A 


1  -  ^2  ^2  |Vomm|2m  2  +  •  •  • 


m  >  1 


(3.18) 


The  terms  retained  on  the  right  side  of  (3.18)  are  independent  of  A,  so  to  this  order 

Cg  =  C. 

We  have  evaluated  V0m  by  numerical  integration -with  1  <  m  <  12,  for  five  values  of 
R.  The  results  for  C(R)/c  from  (3.18)  are 

C(0. 5)  C(0.4)  C(0.3)  C(0.2)  C(0.1) 

=  0.87,  =  0.966,  =  0.991,  =  0.999,  =  1.000. 


(3.19) 


These  results  can  be  compared  with  those  obtained  by  direct  numerical  solution  of  (2.1  )- 
(2.3),  which  are  listed  in  Table  1.  We  see  that  for  small  K  the  agreement  is  good  for  all 
values  of  A,  while  for  small  R  the  agreement  is  good  for  all  values  of  A'.  In  fact  since 
Vmn  =  0(R3),  (3.18)  shows  that  C(R)/c  =  1  —  0(R6),  so  C(R)/c  differs  little  from  unity 
for  R  small.  This  is  in  agreement  with  the  values  in  (3.19)  and  in  Table  1. 

The  matrix  elements  Vmn  given  by  (3.14)  can  be  evaluated  approximately  when  R  is 
small,  which  may  be  useful.  First  we  integrate  by  parts  in  (3.14)  to  get 


Vmn  =  27rj(m 


r « 

-n) 

J-R 


(log  S)e 


2iri(  m  —  n)r 


dx. 


(3.20) 


In  (3.20)  we  write  log  5  =  -/  —  /2/2  —  /3/3  +  . . .  where  /( x)  =  ir(R2  -  x2),  and  obtain 


rR 

Vmn  =  -27ri(m  -  n)  /  f(x)e2*'{m~n)z  dx  +  0(i?5). 
J-R 


(3.2i : 


Integration  yields,  with  amn  =  2n(m  —  n), 


V  = 

v  mn  — 


4t 


- sin(amnf?)  -  Rcos(amnR ) 


+  0(i?5). 


(3.22) 


From  (3.22)  we  see  that  when  |amn.Rj  <  1  then  Vm„  =  amnv  +  0(.R5)  where  v  =  is 

the  volume  of  a  single  sphere. 

The  results  (3.13)  and  (3.15)  are  not  valid  when  any  denominator  in  them  vanishes. 
This  occurs  when  K  —  mr  for  some  integer  n.  Then  ko  =  nr  and  the  general  solution 
of  (3.2)  satisfying  (2.3)  is  po  =  cie,ka 1  +  c2e-,i°x  where  ci  and  c2  are  constants.  To  find 
(fc2)i  we  multiply  (3.3)  by  e±tk°z  and  integrate  from  —1/2  to  1/2.  The  integral  of  the  left 
side  vanishes  yielding  two  equations: 

(k2)ici  -  ikoV’0c2  =  0, 
ik0Vn0ci  +  (k2)iC2  =  0. 

For  C\  and  c2  to  be  non-zero,  the  determinant  of  their  coefficients  must  vanish,  from  which 
we  get  the  two  solutions 

(3.24) 


(3.23) 


(k2)±  =  ±k0\Vn0\. 


Since  S(x)  is  symmetric,  Vno  is  imaginary.  Then  for  (k2) j  =  (fc2)*  (3.23)  yields  ci/c2  =  1 
and  hence  po  =  Pq  =  e,koX(l-f  e,2Tnx).  Similarly  for  (fc2)i  =  ( k 2)f  (3.23)  gives  c2/cj  =  -1 
and  po  =  Pq  e,fc°x(l  -  e~'2*nz). 

Next,  we  calculate  p2.  We  write  pi  as  e'k°zv\(x)  where  iq(x)  has  period  1  and  the 
Fourier  expansion  Ylm  where  the  am  are  constants.  To  find  a;  we  substitute 

Pi  =  e,k°z{Ylm  amelZfTni )  into  (3.3),  multiply  by  e-'(*o+2*-/)  and  integrate  the  resulting 
equation  from  -1/2  to  1/2.  In  this  way  we  obtain 


pf(x)  =  -'k0e 


tkgX 


E 


m^O 


(^mQ  ~F  ^n+m.o) 
[k%  -  (27rm  -f  fc0)2 


(3.25) 


To  calculate  (k2) 2  we  substitute  (3.25)  into  (3.4),  multiply  by  e  ,Ao1  and  integrate 
from  — £/ 2  to  L/ 2.  The  integral  of  the  left  side  vanishes  and  we  obtain  for  Kn  =  rn/L , 

<*’>?  =  j(  E  E  ^ 


m^O,  — n 


Since  Vmo  =  i  T  sin(2xmx)dx  we  have  V£0  =  V-m0.  Also,  for  m  =  —  (/  +  n), 
(n+m)m  ^'m0 m .o  =  —  7— 777^)0 K»+/,o-  Hence  the  second  term  in  (3.26)  vanishes  and 


(3.26)  becomes 


(*•>?- =  E 

4  '  ( n  +  mjm 

m^O.  —  n  '  ' 


We  now  use  the  results  (3.24)  for  (k2)f  and  (3.27)  for  (k2)f  in  (3.1)  for  k 2  with 
kQ  =  K  =  nx.  Denoting  fc  by  kn  we  get 


(*’.)*■=  O'*)2  ±»*iv„.l  +  2  £  £r"i;lv"°l 

4  * — '  ( n  +  m  in 

m/O.-n  V  ’ 


(3.2S) 


The  two  positive  values  k+  an d  determined  by  (3.28)  are  the  two  boundaries  of  the  n-th 

stop  band  corresponding  to  K  =  nx.  the  width  of  this  band  is  A kn  =  k+  —  k~  =  |  Vn0|  -f - 

This  width  is  2irnv  +  •  •  •  when  2 -nR  is  small.  The  values  of  x_1&*  and  of  x-1  A kn  are 
shown  in  Table  3,  based  upon  (3.28).  The  Vm0  were  calculated  by  numerical  integration 
for  m  <  12  for  five  values  of  R  for  the  first  three  bands,  n  =  1, 2,  3.  The  results  agree  well 
with  those  in  Table  2,  obtained  by  numerical  solution  of  the  eigenvalue  problem  (4.1)  with 
IV  +  1  =  100. 


4.  Numerical  solutions 

To  solve  the  eigenvalue  problem  (2. 1-2.3)  numerically  we  discretize  (2.1)  on  {0  <  x  <  L) 
using  a  uniform  grid  Gh ,  where  h  is  the  meshsize.  The  resulting  second  order  accurate 
conservative  discretization  scheme  for  (2. 1-2.3)  is 

~  P? )  -  ~  P.-i)  +  =  0  (4.1) 
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pb-e'KLph,  S_^=S.. 


Ph  „  — **/tL  nh 

_i  =  e  0/-i> 

Here  P/  is  the  solution  at  the  gridpoint  x }  =  jh,  Sa  =  5( ah)  and  A’  +  1  is  the  number  of 
gridpoints  in  Gh .  For  any  given  K  there  are  A'  —  1  solutions  of  the  eigenvalue  problem 
(4.1).  In  our  calculations  we  took  L  —  1  and  A'  —  1  -=  100. 

In  Table  1,  the  numerical  values  of  kh(R)/K  are  given  for  the  first  two  bands,  cor¬ 
responding  to  K/pi  in  the  ranges  0  to  1  and  1  to  2  respectively.  We  have  seen  that  the 
perturbation  results  (3.19)  are  in  good  agreement,  with  them  for  K  small  and  also  for  R 
small. 

In  Table  2,  x-1  [&£ (P)]i  is  given  for  n  =  1,  2,  3  where  [&*]*  axe  the  two  values  of  kh  at 
the  ends  of  the  n-th  stop  band.  For  comparison,  Table  3  shows  the  corresponding  values  of 
k  given  by  (3.28),  obtained  by  second  order  perturbation  analysis,  which  agree  well  with 
the  values  in  Table  2. 

5.  Conclusion 

We  have  presented  two  methods  for  calculating  the  pressure  p(x),  the  phase  velocity  C  and 
the  propagation  constant  K  for  sound  waves  in  a  medium  containing  fixed  rigid  spheres 
arranged  in  a  simple  cubic  lattice.  One  is  a  perturbation  expansion  in  R/L,  where  R  is  the 
radius  of  a  sphere  and  L  is  the  distance  between  centers.  The  other  is  a  direct  numerical 
method.  The  results  of  the  two  methods  agree  well  both  when  R/L  is  small  and  also  when 
KL  is  small. 

The  results  show  that  there  are  pass  bands  and  stop  bands  along  the  axis  of  kL  = 
u-L/c.  The  boundaries  of  the  first  three  stop  bands  are  given  in  Tables  2  and  3  for  five 
values  of  R/L.  The  values  of  C/c  in  the  first  two  pass  bands  are  given  in  Table  1  for 
five  values  of  R/L  and  ten  values  of  KL/ir.  They  show  that  C/c  <  1  in  the  first  pass 

i  .  : 

I 

band,  and  that  C/c  decreases  as  R/L  increases  and  as  KL  increases.  In  the  second  band 
C/c  >  1  except  near  the  upper  edge  of  the  band.  Within  this  band  C/c  still  decreases  as 
KL  increases,  but  it  increases  as  R/L  increases. 
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A 

m 


*~'k(2) 
7T-1  Afc(2 

’r"1*l2) 

TT"1^) 
7r-1  Akir 


7T  lfcJj 


i?  =  0.1 

R  =  0.2 

R  =  0.3 

R  =  0.4 

R  =  0., 

0.9959 

0.9698 

0.9117 

0.827 

0.697 

1.0041 

1.029S 

1.0907 

1.190 

1.366 

0.0082 

0.06 

0.179- 

0.363 

0.669 

1.9927 

1.9652 

1.9637 

2.0099 

2.036 

2.0072 

2.0361 

2.0462 

2.0341 

2.116 

0.0145 

0.0709 

0.4092 

0.0242 

o.os 

2.990 

2.9855 

2.976 

2.9992 

2. 985 

3.0088 

3.0166 

3.0314 

3.025 

3.104 

0.0188 

0.0311 

0.0554 

0.0258 

0.119 

Table  3. 

Boundaries  and  widths  of  the  first  three  stop  bands  calculated  by  second 
order  perturbation  theory,  (3. 28).  Numerical  integration  was  used  for  cal¬ 
culating  the  Vm0,  and  only  terms  with  m  <  12  were  included. 


$ 


